#include <stdio.h>
#include <string.h>
#include <math.h>
#include "userstruct.h"

#pragma region 本算法函数声明_07root
int    Read07EPHFromFile(FILE *hReadFile);
void   clac07root();//计算ECEF数据生成文件
#pragma endregion

#pragma region 外部公共函数声明
int    GetReadCGS2000Para(CGS2000_PAR  *mCGS2000Para);
double Frac(double x);//Fractional part of a number (y=x-[x])
double Modulo(double x, double y);//x mod y
double EccAnom (double M, double e);
double delt_WeekSecT(WeekSecTimeA *_t1,WeekSecTimeA *_t2);//计算两个WeekSecTimeA之间的秒数差
void   WeekSecTimeA_validity(WeekSecTimeA *_tm);/*将时间表达形式规范化*/
#pragma endregion

//<-----------全局变量-----------------
X32          mBrdEPH;          //星历存储变量， 注意：22参数仅用18参数
CGS2000_PAR  mCGS2000Para;     //CGS2000坐标系参数定义
double       FitArcLength=3600;//总拟合弧长3600秒--60分钟
//-----------全局变量----------------->

int main(int argc, char *argv[])
{
	printf("hello c world.\n");
	memset(&mBrdEPH, 0, sizeof(mBrdEPH)); //初始化星历参数初值全部为0

	//1-----读测试根数文件到全局变量中--mBrdEPH
	FILE *hEphFile;
	hEphFile=fopen("./root07ANSI.txt","r");
    if(Read07EPHFromFile(hEphFile))
	{
		printf("open input file error and exit!\n");return -1;
	}
    if(hEphFile)fclose(hEphFile);
    
	//2-----下面计算18根数生成位置文件
    clac07root();
    printf("calc finished.\n");

    return 0;
}

/*
    使用全局变量 mBrdEPH 保存读取的星历文件
    星历文件由提供程序生成
*/
int  Read07EPHFromFile(FILE *hReadFile)
{
    char    ReadChars[300];
    char    TempChars[100];
    char    ReadStr[300];
    int     NoGroup=0;
    int     TempInt=0;
    if (hReadFile == NULL)return 1;
    fseek(hReadFile, SEEK_SET, 0 );//置于文件的开头
    //例子文件7根数含20组数据
    while(!feof(hReadFile) )
    {
        fgets(ReadChars, 300,  hReadFile );
        strcpy(ReadStr,ReadChars);
        //printf("%s\n",ReadChars);  //ReadChars是ANSI编码
		char *subStr1 = "FitArc";
		char *subStr2 = "**";
		char *ch1     = strstr(ReadStr, subStr1);
		char *ch2     = strstr(ReadStr, subStr2);
        if (ch1 != NULL)
        {
            printf("find the id in 7root file.\n");
            sscanf(ReadStr," %s %lf ",TempChars, &FitArcLength);
            FitArcLength *= 60.0;//min to second
        }
        else if (ch2 != NULL)
        {
            printf("ReadStr:%s\n",ReadStr);//** i 869  360004.000000
            sscanf(ReadStr," %s %d %d %lf", TempChars, &TempInt,    //TempChars, &TempInt用于格式化输入参数
                    &mBrdEPH.RNSSbrd_arry[NoGroup].ref_epoch_week,  //周
                    &mBrdEPH.RNSSbrd_arry[NoGroup].ref_epoch_sow);  //星历参考历元，秒，实际上这个数也是toe
            
            //下面是1组7根数的数据
            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].a);

            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].e);

            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].i0);

            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].OMEGA0);

            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].omega);

            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].M0);

            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].idot);

            fgets(ReadChars, 300,  hReadFile);
            strcpy(ReadStr,ReadChars);
            sscanf(ReadStr," %s %lf ", TempChars, &mBrdEPH.RNSSbrd_arry[NoGroup].ref_epoch_sow);

            NoGroup++;
        }
    }

    mBrdEPH.Ngp = NoGroup;//星历组数 <=24
    printf("mBrdEPH.Ngp=%d\n",mBrdEPH.Ngp);//sat12root.rst文件一共有1组12根数数据
    printf("FitArcLength=%f\n",FitArcLength);
    printf("read ok\n");
    return 0;
}

//计算ECEF数据生成文件
void clac07root()
{
    /*时间变量，用于标识需要计算卫星状态信息*/
    WeekSecTimeA       GoalTime;
    char               OutStrChars[200];
    /*用于保存CGS2000的有关参数的定义常量*/
    double  GoalSatPos[3];         //保存目标时刻卫星的位置
    double  CurEphPara[7];     //使用的卫星广播星历，矢量表达
    int   i,nn; //临时变量
    FILE *hRstFile;

    /*时间变量：周＋秒，保存当前星历的参考时间*/
    WeekSecTimeA  CurEphRefTime;

    // Variables
    double  a,e,Parai,Omega,M,M0,n,p,omega;
    double  E,Omegak,tk;
    double  Vk,fik;
    double  uk,rk,ik,xk,yk,tempf;
    double  sinE,cosE,sqrt1_e2;

    GetReadCGS2000Para( &mCGS2000Para );
    if ( (hRstFile = fopen( "ECEF_EPH7Para_POSVEL.rst","wt") ) == NULL )
       return ;

    //临时设置的计算卫星位置的时间起点
    GoalTime.Week = mBrdEPH.RNSSbrd_arry[0].ref_epoch_week;
    GoalTime.Secs = mBrdEPH.RNSSbrd_arry[0].ref_epoch_sow-10*60;

    //如果没有有效的星历，
    if ( mBrdEPH.Ngp <= 0 )
         return ;

    /*当前星历参考时间的周、秒*/
    CurEphRefTime.Week =  mBrdEPH.RNSSbrd_arry[0].ref_epoch_week;
    CurEphRefTime.Secs =  mBrdEPH.RNSSbrd_arry[0].ref_epoch_sow;
    WeekSecTimeA_validity(&CurEphRefTime);//CurEphRefTime.validity();

    i = 0;
    CurEphPara[i++] = mBrdEPH.RNSSbrd_arry[0].a;
    CurEphPara[i++] = mBrdEPH.RNSSbrd_arry[0].e;
    CurEphPara[i++] = mBrdEPH.RNSSbrd_arry[0].i0;
    CurEphPara[i++] = mBrdEPH.RNSSbrd_arry[0].OMEGA0;
    CurEphPara[i++] = mBrdEPH.RNSSbrd_arry[0].omega;
    CurEphPara[i++] = mBrdEPH.RNSSbrd_arry[0].M0;
    CurEphPara[i++] = mBrdEPH.RNSSbrd_arry[0].idot;
    //计算21个点
    for(nn=0; nn<=20; nn++ )   //60秒计算一个点
    { 
        tk = delt_WeekSecT(&GoalTime,&CurEphRefTime);//tk = GoalTime - CurEphRefTime;
        a = CurEphPara[0];
        e = CurEphPara[1] ;
        Parai  = CurEphPara[2]+CurEphPara[6]*tk;
        Omegak = CurEphPara[3];
        omega  = CurEphPara[4];
        n = sqrt(mCGS2000Para.GM/(a*a*a));
        M = CurEphPara[5] + n*tk;
        p     = a*( 1.0-e*e );
        E  = EccAnom(M,e);
        Vk = 2.0*atan( sqrt((1+e)/(1-e))*tan(E/2.0) );  // old from clc
        fik = Vk + omega;
        uk = fik;
        rk = a*(1.0-e*cos(E));
        ik = Parai;
        xk = rk*cos(uk);
        yk = rk*sin(uk);
        Omegak -= mCGS2000Para.omega_earth*tk;
        //计算卫星在地心固定坐标系中的直角坐标
        GoalSatPos[0] = xk*cos( Omegak ) - yk *cos( ik ) * sin( Omegak );
        GoalSatPos[1] = xk*sin( Omegak ) + yk *cos( ik ) * cos( Omegak );
        GoalSatPos[2] = yk*sin( ik );
        //sprintf( OutStrChars,"%12.5E %03d  %7.2lf  %12.3lf  %12.3lf %12.3lf",
        sprintf( OutStrChars,"%09d   %06d   %7.2lf   %12.3lf   %12.3lf   %12.3lf",
                 nn*60, GoalTime.Week, GoalTime.Secs, GoalSatPos[0], GoalSatPos[1],GoalSatPos[2] );
        fprintf(  hRstFile, " %s \n", OutStrChars );
        GoalTime.Secs += 60;  //一分钟算一个结果
        WeekSecTimeA_validity(&GoalTime);//GoalTime.validity();
   }

   if( hRstFile )fclose( hRstFile );

}

